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[ 1 ] Two distinct snowfall events are observed over the region near the Great Lakes 
during 19-23 January 2007 under the intensive measurement campaign of the Canadian 
CloudSat/CALIPSO validation project (C3VP). These events are numerically investigated 
using the Weather Research and Forecasting model coupled with a spectral bin 
microphysics (WRF-SBM) scheme that allows a smooth calculation of riming process 
by predicting the rimed mass fraction on snow aggregates. The fundamental structures of 
the observed two snowfall systems are distinctly characterized by a localized intense 
lake-effect snowstorm in one case and a widely distributed moderate snowfall by the 
synoptic-scale system in another case. Furthermore, the observed microphysical structures 
are distinguished by differences in bulk density of solid-phase particles, which are probably 
linked to the presence or absence of supercooled droplets. The WRF-SBM coupled with 
Goddard Satellite Data Simulator Unit (G-SDSU) has successfully simulated these 
distinctive structures in the three-dimensional weather prediction run with a horizontal 
resolution of 1 km. In particular, riming on snow aggregates by supercooled droplets is 
considered to be of importance in reproducing the specialized microphysical structures in 
the case studies. Additional sensitivity tests for the lake-effect snowstorm case are 
conducted utilizing different planetary boundary layer (PBL) models or the same SBM but 
without the riming process. The PBL process has a large impact on determining the cloud 
microphysical structure of the lake-effect snowstorm as well as the surface precipitation 
pattern, whereas the riming process has little influence on the surface precipitation because 
of the small height of the system. 
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1. Introduction 

[ 2 ] Microphysics of solid-phase clouds and precipitation 
has a very complicated structure compared to liquid-phase, 
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because of a great variety of the particle habits. Remotely 
sensed measurement of ice and snow particles often involves 
large uncertainty due to the various particle densities and 
nonsphericity. Numerical modeling of ice microphysics has 
difficulty in representing their complicated habits of parti- 
cles and various growth mechanisms. Limited in situ mea- 
surements of aircraft or ground-based instruments are useful 
in order to capture the microphysical structure of clouds 
and precipitation and can provide important validation of 
numerical model experiments. 

[ 3 ] Solid-phase particles are empirically classified into 
distinct categories such as cloud ice, snow, graupel and/or 
hail in a numerical modeling framework of cloud micro- 
physics [e.g., Lin et al., 1983], according to the difference in 
terminal fall velocity. These models involve spontaneous 
autoconversion from one category to another. A few 
microphysical models, however, have introduced advanced 
approaches to cloud microphysical representations using a 
gradation of distinct categories: A multicomponent (water 
mass, solute mass and ice shape factor) bin model was 
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developed using a hybrid bin method [Chen and Lamb, 
1994]; ice volume was recently introduced in the method- 
ology to represent continuous changes in the bulk density of 
particles [Misumi et ah, 2010]. A spectral habit ice predic- 
tion system (SHIPS) that can retrain the history of particle 
growth by predicting particle property variables (PPVs) was 
developed [Hashino and Tripoli, 2007]. A conceptual model 
considering mass-dimension and area-dimension relation- 
ships as a continuous function of particle size and rimed 
mass fraction was integrated into a bulk microphysical 
model [Morrison and Grabowski, 2008] and into a bin 
microphysical model [Morrison and Grabowski, 2010], 
Snow and graupel were included in the precipitation ice 
(PI) category using a varying rime intensity parameter in the 
traditional bulk microphysical parameterization [Lin and 
Colle, 2011; Lin et al., 2011], Prediction of rimed mass in 
each of snow aggregates bins was introduced into a typical 
1 moment spectral bin microphysics (SBM) to allow a 
smooth treatment of the transformation from snow to 
graupel or hail [Khain et al., 2011, 2012], 

[ 4 ] The scientific goal of this study is to research the fol- 
lowing points: Although these advanced approaches have 
been developed, there still remains a lack of full under- 
standing of how ice cloud microphysics is characterized in 
the real atmosphere and how the numerical simulation can 
be optimized using the advanced microphysical approaches 
noted above. A comparison between numerical simulations 
and observational data is of considerable importance to 
confirm the reproducibility of the simulation and to investi- 
gate factors that lead to differences in the microphysical 
structure. 

[5] This paper focuses on two particular snowfall events 
on 19-23 January 2007 observed during the field campaign 
of the Canadian CloudSat/Cloud-Aerosol Lidar and Infrared 
Pathfinder Satellite Observation (CALIPSO) Validation 
Project (C3VP). The campaign was conducted originally 
for the purpose of validating CloudSat/CALIPSO retrieval 
algorithm using both in situ and remotely sensed obser- 
vations. The concentrated field measurements were carried 
out in south central Ontario, Canada from 1 October 2006 
to 31 March 2007. The Center for Atmospheric Research 
Experiments (CARE) is the main site of the ground-based in 
situ and remotely sensed measurements including a lidar, 
radars, microwave radiometers, a spectrometer, precipitation 
gauges, and disdrometers. In addition, flight measurements, 
containing a suite of in situ and remotely sensed observa- 
tions, were conducted around the CARE site. 

[6] Shi et al. [2010] conducted numerical simulations of 
the two snowfall events using the Weather and Research 
Forecasting (WRF) model with a newly implemented 
Goddard microphysics scheme (1 moment bulk for 2 water 
and 3 ice classes). Their simulations successfully reproduced 
the radar reflectivity distributions of the snowfalls observed 
by the King City C-band radar measurements. The vertical 
cross section of the cloud structures also agreed well with 
corresponding CloudSat observations. The Shi et al. study 
mostly focused on the cloud macrophysical structures of the 
snowfall events using various observational data; however, 
the cloud microphysics was not explored in detail and 
remains to be fully investigated. 

[7] This study utilizes a newly developed WRF coupled 
with spectral bin microphysics (WRF-SBM) that is identical 


to the microphysical part of the Hebrew University Cloud 
Model (HUCM) [Khain et al, 2011, 2012], The WRF- 
SBM is advantageous for studies involving comparison with 
in situ microphysical measurements, because it allows an 
explicit calculation of cloud particle size distributions 
(PSD), which are important to determine the microphysical 
properties of clouds and precipitation. In addition, the 
updated SBM can diagnoses changes of bulk density of 
snow aggregates by the explicit prediction of rimed mass 
fraction on snow. This function is of importance to investi- 
gate impacts of riming process on ice cloud microphysics. 

[8] The Goddard satellite data simulator unit (G-SDSU) 
[Matsui et al., 2010] has been introduced to enable direct 
comparison between the WRF-SBM simulation and in situ 
or remotely sensed measurements by modulating the simu- 
lation output to be compatible with the observations. Using 
the G-SDSU, interactive analysis through the three compo- 
nents (remotely sensed, in situ measurements and model 
simulation) is possible to facilitate understanding of the 
complexities of solid-phase cloud microphysics. 

[ 9 ] In this paper, the methodology including the outline of 
the simulation is described in section 2. The results of 
numerical experiments and comparison with measurements 
are discussed in section 3 . The summary and conclusions are 
given in section 4. 

2. Methodology 

2.1. Description of the WRF-SBM 

[ 10 ] The WRF model is a regional numerical weather 
prediction (NWP) system principally developed by the 
National Center for Atmospheric Research (NCAR) in col- 
laboration with several research institutions in U.S. The 
WRF has two types of dynamical solvers: One is the 
Eulerian mass solver referred to as the Advanced Research 
WRF (ARW). The other is the Nonhydrostatic Mesoscale 
Model (NMM) developed by the National Centers for 
Environmental Prediction (NCEP). The ARW version 3.1.1 
is employed in this study with the following physical 
options: The eta surface layer scheme based on similarity 
theory [Monin and Obukhov, 1954] is used to calculate the 
surface boundary layer dynamics, and Mellar- Yamada- 
Janjic Level 2.5 turbulent closure model [e.g., Janjic, 1990] 
is employed to calculate the planetary boundary layer (PBL) 
process. The Noah land surface model with a four-layer soil 
is applied to provide heat and moisture fluxes over land 
surface. The updated Goddard radiation package for both 
longwave and shortwave radiations is used to compute 
atmospheric heating and cooling rate profiles and surface 
energy budget. 

[ 11 ] The WRF ARW version 3.1.1 was coupled with the 
SBM part of the HUCM [Khain et al., 2011]. This version of 
SBM is the full package that is different from the Fast-SBM 
described by Lynn et al. [2005a, 2005b] and Khain et al. 
[2009, 2010], Cloud hydrometeors are categorized into 
one-water and six-ice classes, that is, water droplets, ice 
crystals (plate, column, dendrite), snow aggregates, graupel, 
and hail. The discrete PSDs of the hydrometeor classes are 
represented on a grid containing 43 doubling mass bins 
covering particles mass sizes in a range of 3.35 x 10 -11 
g < mass < 1.47 x 10 2 g (2 /jm < radius <32.8 mm in terms 
of the radii of droplets or melted ice). The relationships 
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(a) 



radius of mass equivalent sphere with the bulk density [urn] 


(b) 



Radius of mass equivalent sphere with the bulk density [urn] 

Figure 1 . Relationships (a) between the particle bulk den- 
sities and radii of mass equivalent sphere with the bulk den- 
sity and (b) between the terminal fall velocity at 1000 hPa 
and radii of mass equivalent sphere with the bulk density 
assumed in the spectral bin microphysical model. 

between the particle bulk densities and radii of mass 
equivalent sphere with the bulk density assumed in the 
SBM are shown in Figure la. Snow aggregates, graupel, 
and hail are assumed to be fluffy spheres when calculating 
their microphysical processes [Khain and Sednev, 1995], 
The supplementary mass size distributions representing a 
liquid part on melting ice particles and a rimed part on snow 
aggregates are also added into the sets of 43 bins. The size 
distribution of condensation nuclei (CN) is discretized into a 
mass grid containing 13 bins with a radius range from 10~ 3 
/jm to 1 n m. The use of a smaller number of CN bins than 
that of the original HUCM SBM is to improve the effi- 
ciency of computation [Iguchi et al., 2008]. 

[ 12 ] The SBM calculates nucleation of droplets and ice 
crystals, including condensation and deposition growth, 
evaporation, sublimation, droplet freezing, riming, melting, 
shedding, coalescence growth, and breakup of the categorized 


hydrometeor particles. The nucleation of droplets is calcu- 
lated on the basis of the Kohler equation using grid-scale 
supersaturation with respect to water. All CN lager than 
critical radius determined by the supersaturation value is 
assumed to be converted to cloud droplets of the radius that 
is determined by the equation. 

[ 13 ] The formation processes of ice crystals are divided 
into homogenous and heterogeneous nucleation and ice 
multiplication. The heterogeneous nucleation is further 
classified into four types: deposition, condensation-freezing, 
immersion-freezing and contact-freezing nucleation [Rogers 
and Yau, 1989]. The SBM explicitly calculates ice multi- 
plication [Hallett and Mossop, 1974] and heterogeneous 
nucleation, except for contact-freezing nucleation. Contact- 
freezing nucleation is neglected in the present model 
because of their lower production efficiencies [ Khain et al., 
2000]. A number concentration of ice nuclei (IN) is not 
treated as an explicit prognostic variable, but its function is 
implicitly parameterized. The equation of Meyers et al. 
[1992] for deposition and condensation-freezing nucleation 
is used to calculate a newly generated number of primary ice 
crystals at each time step; the diagnostic equation is formu- 
lated into the time integration using a semi-Lagrangian 
approach [ Khain et al., 2000]. The type of generated ice 
crystal is dependent on the environmental temperature 
[Takahashi et al., 1991]. Immersion-freezing nucleation is 
calculated using two types of parameterizations: When the 
ambient temperature is higher than 235.15 K, the probability 
of freezing for supercooled droplets is calculated on the 
basis of the parameterization of Ovtchinnikov and Kogan 
[2000], Otherwise, the parameterization by Bigg [1953] is 
used to provide the probability of freezing of supercooled 
droplets. Homogenous nucleation, that is, homogenous droplet 
freezing, is currently included in the Bigg’s parameterization. 

[ 14 ] Condensation, deposition, evaporation, and sublima- 
tion are calculated using the solution proposed by Khain 
et al. [2008]. A flux method of Bott [1998] is employed 
for solving the stochastic collection equation to determine 
the coalescence growth. This scheme is diffusive and hence 
the particle number concentration and the mixing ratio are 
potentially underpredicted [ Connolly et al., 2012]; the 
problem should be improved in future work. The collision 
efficiencies are given from the studies of Ji and Wang 
[1990], Wang and Ji [1997], Pruppacher and Klett [1997], 
Pinsky et al. [2001], Khain et al. [2001], and Khain et al. 
[2004]. Gravitational sedimentation is calculated using a 
box-Lagrangian drop scheme [Kato, 1995]. The terminal fall 
velocities are given from the studies of Beard [1976], 
Pruppacher and Klett [1997], and Khain et al. [2001]. 
Gradual melting of snow aggregates, graupel and hail and 
shedding of water from melting graupel and hail are calcu- 
lated using the supplementary size distribution of liquid 
mass inside the particles; the methodology of the calculation 
follows Rasmussen and Heymsfield [1987] and Phillips etal. 
[2007], 

[ 15 ] Note that the SBM can explicitly calculate changes of 
bulk density and resultant bulk radius of snow aggregates 
using the supplementary mass size distribution of rimed part 
in each mass bin at every time step. If the bulk density of 
aggregates with bulk radius larger than 550 pm exceeds 
0.2 g/cm 3 , the content of the bin is moved into the same 
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mass bin of graupel. In addition, if the rimed mass fraction 
of aggregates with bulk radius larger than 550 pm is over 
90%, the similar migration to graupel then hail occurs. The 
terminal fall velocities and the collision kernels of snow 
aggregates are corrected at any time according to changes of 
the bulk density and radius. The terminal velocities of rimed 
aggregates are calculated by a liner interpolation between the 
values of nonrimed snow aggregate and hail. The collision 
kernels are recalculated using the modified radii as well as 
terminal fall velocities. 

2.2. Setup of Numerical Experiments 

[16] Two particular snowfall events from 19 to 23 January 
2007 are the focus of this study; the first and second events 
are referred as to lake-effect snowstorm and synoptic 
snowfall, respectively, according to Shi et al. [2010], Two 
36 h individual integrations were conducted from 12:00 UTC 
on 19 January to 00:00 UTC on 21 January, and 12:00 UTC 
on 21 January to 00:00 UTC on 23 January. Three simula- 
tion domains were constructed with horizontal grid intervals 
of 9, 3, and 1 km and corresponding grid components of 
301 x 241, 430 x 412, and 457 x 457 for the dOl, d02, and 
d03 domains following those in the work of Shi et al. [2010, 
Figure 4]; the inner domain d03 is centered at 44.0°N, 80.7° W. 
The vertical domain, extending to a height of approximately 
20 km was divided into 60 layers with intervals increasing 
with altitude. 

[ 17 ] First, the domain dOl is individually simulated using 
the initial and lateral boundary conditions from the analysis 
data sets of the North America Mesoscale model forecast 
(NAM) 218 AWIPS grids; the data sets are distributed by 
NCEP with 12 km horizontal grid intervals, 40 vertical levels 
and 4 samples per day. The one-moment bulk micro- 
physics of Goddard Cumulus Ensemble model (GCE) \Tao 
and Simpson , 1993; Tao et al., 2003; Lang et al., 2007] is 
employed in this domain, together with the Grell-Devenyi 
ensemble cumulus parameterization [Grell and Devenyi, 
2002], The 3ICE-graupel option [Tao et al., 2003] includ- 
ing cloud water, rain, cloud ice, snow, and graupel is 
selected. 

[is] Second, the domains d02 and d03 are simulated at the 
same time using online two-way nesting; all prognostic 
variables are common between domains d02 and d03 in the 
online two-way nesting configuration. An offline one-way 
nesting is applied to connect the domain d02 to the parent 
domain dOl . No hydrometeor advection is considered between 
domains dOl and d02 in the one-way nesting configuration. 
The initial conditions of d02 and d03 are given from the 
analysis data sets of the NAM 218 AWIPS grids. The fully 
packaged SBM is employed in both domains d02 and d03 
without a subgrid cumulus parameterization. The initial CN 
concentration field is uniform in the layer under a height of 
2 km and decreases exponentially with a scale height of 
2 km above 2 km. The initial CN size distribution is set 
using an empirical relationship between the number con- 
centration of cloud condensation nuclei (CCN) and super- 
saturation with respect to water [Khain et al., 2000, 
equations 2.4 and 2.6]: 

ALpi.5^1.5 ^(TU) W , (1) 


where No and k are the parameters that determine the con- 
centration and the soluble fraction of CN in the assumed air 
mass. No = 1500 cm~ 3 and k = 0.5 were set in this experi- 
ment; these values were selected by consideration of the 
local environment of continental area near the lakes [e.g., 
Pruppacher and Klett, 1997; Khain et al., 2010], A and B 
are the coefficients in the Kohler equation [e.g., Rogers and 
Yau, 1989] considering the Kelvin effect and Raoult’s Law, 
respectively. Zero horizontal gradient was set for the inflow 
case in the lateral boundary region of domain d02 to prevent 
cleanup of CN by incoming advection. 

2.3. Goddard Satellite Data Simulator Unit 

[ 19 ] The G-SDSU is an off-line module unit to convert 
output of atmospheric model simulation into the variables 
that correspond to primary products of various measure- 
ments. The output can be used for the evaluation of simula- 
tions through direct comparison with low-level measurement 
products [ Matsui et al., 2009; Masunaga et al., 2010]. This 
approach for comparison between simulations and measure- 
ments has the following advantage: Retrieval algorithms 
generally include some assumptions to make geophysical 
parameters from the raw products. The assumptions are often 
in disagreement with those used in the atmospheric model. 
This inconsistency can make the comparison unreliable, 
when the geophysical parameters derived from the retrieval 
and simulation are compared. However, the direct compari- 
son using the simulator unit can bypass the problem, because 
the assumptions in the simulator unit are consistent with 
those made in the atmospheric model. 

[ 20 ] The G-SDSU is run off-line and coupled with the 
WRF-SBM directly to make use of simulated PSDs of 
hydrometeors and rimed mass size distribution of snow 
aggregates. Although many kinds of measurement activities 
were conducted in the C3VP field campaign, this study has 
mostly focused on the observational data sets of King City 
C-band radar, NRC Convair 580 aircraft and ground-based 
in situ instruments at CARE site. The methods to calculate 
the primary products from the WRF-SBM simulation are 
explained in section 3, together with an analysis of the 
comparative results. 

3. Results 

3.1. Analysis of the Macro Profiles Using C-Band 
Radar Data 

[ 21 ] The King City 5.5 GHz C-band dual-polarization 
radar is located at 43.96°N, 79.57°W, 30 km from the CARE 
site (44.23°N, 79.78°W). During C3VP campaign, the radar 
performed it’s normal operation scan strategy, collecting 
multiple tilt data every 10 min. Special single tilt scans 
tailored to C3VP requirements were inserted into the 
sequence; however, for the purpose of this study only the 
multiple tilt data were used. The details of the radar char- 
acteristics can be found in the work of Hudak et al. [2006] 
and Huang et al. [2010]. Subsequent to quality control 
procedures, the scanning radar reflectivity was interpolated 
to a three-dimensional Cartesian coordinate grid with a 
200 km domain extending north-south and east-west. The 
grid resolution was 1 km and 500 m in the horizontal and 
vertical, respectively. This study used only gridded copo- 
larized radar reflectivity. The reflectivity data is processed 
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to compute the maximum reflectivity (Z max ) and 0 dB echo- 
top height (H et ). 

[ 22 ] The corresponding product of the C-band radar is 
calculated using the output of the WRF-SBM simulation 
through G-SDSU. The methodology to calculate radar 
reflectivity in the G-SDSU is based on that of Masunaga 
and Kummerow [2005]. The backscattering cross section of 
particles is calculated on the basis of the full solution of 
Mie-based routine with dialectic constant through the 
Maxwell-Garnet assumption (ice inclusion with air matrix) 
and assuming the particles to be a spherical body with the 
particle density and size assumed in the SBM for each mass 
bin. The spherical assumption is highly appropriate to the 
calculation of the microwave frequency range (C-band) in 
Rayleigh regime [Liu, 2008]. 

[ 23 ] Figures 2 and 3 show the spatial distributions of Z max 
and H et for the lake-effect and synoptic events in the form of 
snapshots at 00Z, 06Z and 12Z on 20 and 22 January. The 
observed Z max patterns of the both events (Figures 2a 
and 2c) correspond reasonably well to those calculated 
from the WRF-SBM simulations through the G-SDSU radar 
simulator (Figures 2b and 2d). In the lake-effect case, an 
isolated intense snowband was triggered by the passage of a 
cold front over Lake Huron. Cold dry air blowing across the 
lake accentuated the heat and moisture fluxes from the rel- 
atively warm water surface of the lake, and formed the 
snowband on the leeward shore. At 00Z on 20 January, the 
formation of the lake-effect snowbands with H et of approx- 
imately 2 km was observed in both the observational data 
and simulation in the vicinity of Geogian Bay of Lake 
Huron. At 06Z, both the observation and simulation show 
two snowbands extending northwest-southeast, though the 
locations of simulated bands are mismatched to the radar 
observations. At 12Z, one snowband in the vicinity of 
Georgian Bay remained. The magnitude of Z max is approx- 
imately 25 dB in both the simulation and observations and is 
distributed almost uniformly from 00Z to 12Z. 

[ 24 ] In contrast to the lake-effect case (20 January), the 
precipitation system on 22 January was synoptically 
driven with widespread moderate snowfall caused by the 
passage of synoptic low-pressure system. At 00Z and 06Z 
on 22 January, a relatively homogenous Z max distribution 
covered almost the entire domain of the snapshot. At 12Z, 
the observation showed an end of precipitation as the sys- 
tem moved eastward, whereas the simulation indicated the 
persistence of weak radar echo at this time. The Z max in the 
synoptic event ranges from approximately 15 dB to 25 dB 
and H et is approximately 5 km with the exception of 12Z, 
where the H et is notably lower as the system decayed. Note 
that the circular pattern with low H et around the center of 
plots at 00Z and 06Z in Figure 3c is an artifact of the 
interpolation process with the limited vertical elevation 
angle of the ground radar. 

[ 25 ] The comparison of the snapshots in Figures 2 and 3 
shows some forecast errors for the location and timing of 
both systems. This kind of error is often controlled by initial 
and boundary conditions rather than the model physics. A 
quantitative analysis that does not depend on location and 
timing can be useful to compare the overall structures of the 
observed and simulated precipitation systems. Figure 4 
shows a statistical comparison in the form of joint altitude- 
reflectivity contoured frequency diagrams, illustrating the 


overall probability density for both the altitude and reflec- 
tivity. The data accumulated from 22Z on 19 January to 12Z 
on 20 January for the lake-effect case and from 00Z to 12Z 
on 22 January for the synoptic case over the domain of 
200 km square centered at the radar is sampled to make the 
statistical plots. The areas with radii of 15 km in the lake- 
effect case and 40 km in the synoptic case around the radar 
are not included to remove the effect of the limited vertical 
sampling of the radar and the subsequent impact on the 
Cartesian interpolation. 

[ 26 ] The simulation plots share common characteristics 
with those of the radar observations and show the difference 
in the profiles between the two events (Figure 4). The 
probability distributions are limited in the vertical to alti- 
tudes below 3 km and encompass a wide range of radar 
reflectivity (30 dB) in the lake-effect case. In contrast, the 
corresponding distributions extend over 4 km in depth and 
are slightly narrower with maximum reflectivities less than 
30 dB in the synoptic case. Together with the analysis of 
Figures 2 and 3, the lake-effect event is characterized as a 
localized small height and intense snowstorm, derived 
principal from lake-effect processes, whereas the synoptic 
event is characterized as a widely distributed, relatively 
high-up with moderate snowfall. 

[ 27 ] In the lake-effect case, the simulated probability dis- 
tribution is biased to the lower values relative to the radar 
observations in terms of both altitude and reflectivity. This is 
considered to be due to limitations of both the measurements 
and simulation. The radar coverage area is generally at an 
elevation of several hundred meters, and the Blue Mountains 
(44.5°N, 80.4°W) with the maximum terrain height of 
approximately 550 m is located near the southern side of 
Geogian Bay (approximately 90 km from the King radar). 
Because the lowest elevation angle utilized in the King City 
radar data is 0.3°, the height of the radar beam is at an 
approximate elevation of 1 .2 km above ground level at this 
range. Therefore, there is some uncertainly in measuring 
and interpolating reflectivity at the lowest heights shown 
in Figure 4. Yet, the simulation still underestimates the 
height of the system as illustrated in Figures 3a and 3b. The 
observed echo-top height is approximately 1 km taller than 
that of the simulation. The similar shallow bias of the sim- 
ulated snowband is also reported by Shi et al. [2010]. 

[ 28 ] In the synoptic case, the simulated distribution is in 
overall better agreement with the observed pattern, though 
the simulated distribution is slightly biased to lower reflec- 
tivities and higher altitudes compared to the radar. We 
speculate that these small biases are due to forecast errors 
that are caused by imperfect initial and boundary conditions 
and not removed in the quantitative analysis because of the 
limitation of sampling. The magnitude of probability density 
is underpredicted, because the distribution profile of the 
simulation is wider in the vertical compared to the radar 
observations. In contrast, the magnitude of the simulated 
probability density is overpredicted in the lake-effect case 
because the simulated distribution is narrower in the vertical 
compared to the observation distribution. 

[ 29 ] Observations at the CARE site showed a distinct 
difference in the surface precipitation amount between the 
two events. Approximately 12.3 mm and 2.4 mm of liquid 
water equivalent (LWE) 24 h accumulated surface precipi- 
tation were measured on 20 (lake-effect case) and 22 Jan. 
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a) Lake-effect, measurement 

00Z20JAN 06Z20JAN 12Z20JAN 




00Z22JAN 


c) Synoptic, measurement 

06Z22JAN 




d) Synoptic, simulation 



-5 0 5 10 15 20 25 30 35 40 

Figure 2. Snapshots of horizontal distributions of vertically maximum C-band radar reflectivity (dB) 
(a, c) derived from King City radar measurement and (b, d) simulated by Weather Research and 
Forecasting model coupled with a spectral bin microphysics (WRF-SBM). 
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d) Synoptic, simulation 
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Figure 4. Normalized contoured frequency with altitude and reflectivity diagram derived from (a, c) King 
City radar measurement and (b, d) WRF-SBM simulation. 


(synoptic case), respectively, by the double-fenced inter- 
national reference (DF1R) gauge at the CARE site [ Huang 
et al., 2010]. The corresponding values in the WRF-SBM 
simulations are 3.39 mm for the lake-effect event and 
2.16 mm for the synoptic event. This significant under- 
prediction of surface precipitation rate in the lake-effect case 
is considered to be mostly due to a forecast error of the 
location of intense snowband. The simulated intense snow- 
fall region is biased toward the western side, where the Blue 
Mountains are located. In contrast, the observed snowband 
occurred over the CARE site; 06Z snapshots of Figures 2a 
and 2b illustrate the difference in observed and simulated 
snowband location. The maximum value in the simulated 
snowband region is more than 15 mm (Figure 9a), which is 
in good agreement with the surface precipitation amount 
observed at the CARE site. 

3.2. Analysis of the Microphysical Structures Using 
In Situ Observational Data 

[ 30 ] In this subsection, in situ aircraft and ground-based 
measurements are focused and compared with the simulated 
results in order to investigate the microphysical structures of 
the two snowfall systems. 


3.2.1. Aircraft Measurement 

[ 31 ] The NRC Convair 580 aircraft, equipped with a series 
of in situ and remote sensor, was employed to conduct the 
flight mission along CloudSat orbit between Ottawa and the 
CARE site during the measurement campaign. The two 
flight legs during 19-23 January 2007 fell predominantly 
within a domain with a radius of ±0.5° latitude and longi- 
tude centered at 79.5°W, 49.5°N. The first and second 
flights continued approximately from 22Z on 19 to 02Z on 
20 January (the lake-effect case) and from 00Z to 06Z on 
22 January (the synoptic case), respectively. Data sampled 
by the instruments was recorded every 5 s. 

[ 32 ] Cloud particle size distribution was measured by the 
Particle Measurement System’s (PMS) two-dimensional 
cloud (2D-C) and two-dimensional precipitation (2D-P) 
imaging probes. Sampled particles were divided into 33 size 
bins with midpoint diameters ranging from 30 pm to 2.5 cm 
(A. J. Heymsfield et al.. Ice cloud particle size distributions 
and pressure-dependent terminal velocities from in situ 
observations at temperatures from 0 to — 86°C, submitted to 
Journal of the Atmospheric Sciences, 2012). Note that the 
sampling does not separate particles into particular hydro- 
meteor categories. Concentrations in bins below 100 pm are 
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not reliable and are not used to estimate the ice water content 
from the particle size distributions. Total condensate water 
content (TWC) was measured by the DMT cloud probe and 
counterflow virtual impactor (CSI). The CS1 probe has 
approximately 11% uncertainty for TWC of 0.2 g/m 3 and 
larger uncertainty for smaller TWC, when operating in its 
normal configuration [Twohy et al., 1997]. It cannot be 
reliably detect TWC below 0.01 g/m 3 ; however, the CSI 
probe has a shorter region for condensate evaporation which 
adds to the uncertainty in a yet uncharacterized way. 

[ 33 ] The bulk effective particle radius (PMS r e ) is defined 
for the PMS 33 bins as the following form: 

33 

X>m 

PMS r e = , (2) 

J2n 2 Ni 

i = 3 

where r, is a half of the maximum diameter of /,/, bins and TV,- is 
the particle number concentrations. By combining CSI mea- 
surements, the bulk effective particle density (CSI-PMS p e ) is 
formulated in the following form [Heymsfield et al., 2004] 
using TWC estimated by the CSI and particle volume con- 
centration assuming that all PMS particles are spherical: 



[ 34 ] These parameters correspond to representatives of 
particle size and density measured by the airborne instru- 
ments. No liquid water is included in the calculation of these 
parameters. Note that CSI underestimates the water content 
calculated by PMS 2D PSD by at most order of 50%, so 
that there is corresponding possibility of underestimation 
of CSI-PMS p e . 

[ 35 ] The corresponding PMS r e and CSI-PMS p e are cal- 
culated from output of the WRF-SBM simulation on the 
basis of the following approach for comparison with the 
aircraft measurement: The simulation output including 
hydrometeors PSDs is recorded per hour (model time). At 
grid points in the output, each 43 binned PSD of the 6 ice 
hydrometeor categories is distributed over the 33 bins 
identical to the binned size ranges of the aircraft PMS pro- 
gram, and then one 33 binned PSD is recomposed by total- 
ing. This is an imitation of the actual PMS sampling to 
clouds. Note that the two-dimensional video disdrometer 
(2DVD) correlation is used to take into account realistic 
irregular shapes of snow aggregates in the distribution cal- 
culation; detail of the correction is summarized in Appendix 
A. The mass-size relationships assumed in the SBM [e.g., 
Khain and Sednev, 1995] are used for the calculation of the 
other hydrometeor categories. Then, PMS r e (equation (2)) 
and the denominator of equation (3) can be calculated using 
the recomposed 33 binned PSD. TWC estimated by the CSI 
(the numerator of equation (3)) is assumed to be identical to 
the total mass concentration of 43 binned PSDs of the 6 ice 
hydrometeor categories and subsequently CSI-PMS p e can 
be calculated using equation (3). The two parameters are 
finally sampled from the model grid points only in the 
vicinity of the cross sections, during the aircraft sampling 


periods. The two flight legs of the lake-effect and synoptic 
cases are roughly approximated to be the cross sections of 
(44.2°N, 80.0°W)-(44.8°N, 79.8°W) and (43.9°N, 79.6°W)- 
(44.6°N, 79.5°W), respectively. Although an effect of the 
instrument noise is not considered, this approach allows a 
proper comparison of microphysical structure between the 
simulation and measurement [Matsui et al., 2009]. 

[36] Figures 5a and 5c show scatter diagrams between 
PMS r e and CSI-PMS p e derived from the WRF-SBM 
simulations and the aircraft measurements. This plot indi- 
cates a bulk size-dependent spectrum of particle density. 
Observed r e -p e relationships (triangle marks) show that p e 
decreases with an increase in r e , generally. This means that 
the maximum dimension to mass tends to be increased 
when ice crystals are involved in aggregation. A comparison 
between the observed relationships of the two cases shows 
that the synoptic case tends to have smaller density with less 
variability as a function of size. This result has led to the 
speculation that the microphysics structure of the synoptic 
case is dominated by the nonrimed elongated snow aggre- 
gates with relatively small density, while the lake-effect case 
tends to have rimed snow aggregates with large density. 

[ 37 ] The corresponding plots of the simulation (square 
marks) show that r e -p e relationships of the simulation and 
measurement have a common point that is distinct between 
the two cases. A group of relatively large density around r e 
of 700 pm is simulated and observed in the lake-effect case, 
as compared with the plots of the synoptic case. This dif- 
ference can be explained by density change through pres- 
ence or absence of riming on snow aggregates by collision 
with supercooled droplets. In the lake-effect case, the micro- 
physical structure of the intense snowstorm is characterized 
by presence of supercooled water and resultant increase in 
particle density through riming process. The simulated r e -p e 
relationship is not a logarithmic line form in which p e 
decreases with an increase in r e . The scatter of plots spreads 
out around r e of 550 pm to p e of approximately 0.4 g/cm 3 , 
which is equal to the bulk density of graupel assumed in the 
SBM (Figure la). This pattern seems to represent a transfer 
from nonrimed snow aggregates to graupel via rimed snow 
aggregates. This result is consistent with the fact that snow 
aggregates and graupel are dominant on the corresponding 
model grid points for the airborne comparison (Figure 6b). 
The mechanism characterizing the simulated r t .-p c relation- 
ship is consistent with the speculation that the observed 
clouds include rimed snow aggregates with a relatively large 
density, though the simulation cannot reproduce a part of 
large r e side over approximately 1000 pm. Figure 6b also 
makes sure of the presence of supercooled water that can be a 
source of riming by collision with snow aggregates in the 
simulation. All liquid droplets are supercooled because 
atmospheric temperature is less than 0°C in all layers in this 
case. The microphysical structure with supercooled water 
and resultant riming can be implicated in the general weather 
condition by the following: The relatively cold and strong 
wind on the lake surface [Shi et al., 2010, Figure 2] causes a 
large upward moisture flux from the lake; the value of 
upward moisture flux is approximately twice as much as in 
the synoptic case on average. Saturation with respect to both 
water and ice is accomplished in the boundary layer clouds of 
low saturation pressure because of the low air temperature. 
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Figure 5. Scatter diagrams (a, c) between airborne instrument-based bulk density and bulk effective 
radius sampled for the water contents larger than 0. 1 g/m 3 and (b, d) between ground instrument-based 
bulk density and bulk effective radius. Triangle and square marks denote the values derived from the mea- 
surement and the simulation, respectively. 


[38] In contrast, the calm snowfall in the synoptic case is 
characterized by nonrimed snow aggregates and absence of 
supercooled water. The simulated r e -p e relationship is loga- 
rithmically linear and p e uniquely decreases with an increase 
in r e , following the corresponding relationship of nonrimed 
snow aggregates assumed in the SBM (Figure la). This result 
is consistent with the fact that snow aggregates are signifi- 
cantly dominant on the sampled grid points (Figure 6d). 
Rimed snow aggregates and graupel with large particle den- 
sities are not generated in this case, because there is no 
supercooled water. The simulated pattern of the r e -p e plots is 
overall agreement with the observed one with exception of 
slight underestimation of r e . The feature of the simulated 
result does not contradict the speculation that the observed 


clouds are occupied by nonrimed snow. In the synoptic case, 
only ice saturation is accomplished in the simulated clouds. 
Moisture supply from the lake has an insignificant effect on 
the snowfall system, because small moisture supply by the 
weak surface wind hardly reach the synoptic system clouds 
aloft. 

[ 39 ] Furthermore, the simulated presence or absence of 
supercooled water can be directly verified in the same air- 
craft measurements. Presence of supercooled droplets is 
estimated using voltage signal measured by the Rosemount 
Icing (RICE) probe [Heymsfield and Miloshevich, 1989]. 
The RICE measurement shows that the anomalous voltage 
signals larger than 2.2 V were observed in the lake-effect 
cases around the altitude ranging from 1 km to 3 km, which 
indicates presence of supercooled droplets (Figure 6a). This 
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Figure 6. (a, c) Vertical distribution of Rosemount Icing (RICE) probe voltage signals sampled by air- 
craft measurements, (b, d) Vertical distribution of the dominant hydrometeor types in the WRF-SBM 
simulation for the corresponding grid points of the airborne sampling. 


result justifies the presence of supercooled water in the 
simulation of the lake-effect case (Figure 6b). In contrast, 
there are no abnormal RICE signals for supercooled water in 
the synoptic case (Figure 6c), so that the absence of super- 
cooled water in the simulation of the case is also justified 
(Figure 6d). 

[ 40 ] In addition, the aircraft measurements may provide a 
sign of graupel together with supercooled water. Figure 7 
illustrates a correlation between RICE voltage signals and 
particle area ratio with a diameter of approximately 800 p m. 
A high value close to 1 .0 of area ratio means that the particle 
has a spherical shape. Area ratios of high magnitudes 
involved in anomalous RICE signals in the lake-effect case 


demonstrate a transition from snow aggregates to spherical 
graupel as a result of riming. In contrast, the parameters of 
slightly lower magnitudes in the synoptic case can be inter- 
preted as continuous existence of nonspherical snow aggre- 
gates without riming. 

[ 41 ] The WRF-SBM simulation cannot reproduce observed 
low p e near 0.01 g/cm 3 . This is due to the relationship 
between density and size of nonrimed snow aggregates 
assumed in the SBM (Figure la). Another highlighted dif- 
ference between the simulation and measurement is under- 
estimation of r e in the lake-effect case. This point is discussed 
with additional simulations in the next sub section 3.3. 
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Figure 7. Time series of ice water contents (IWC) calculated from the Particle Measurement System’s 
(PMS) two-dimensional cloud (2D-C) and two-dimensional precipitation (2D-P) imaging probes and 
counterflow virtual impactor (CSI), RICE voltage signal, and area ratio in the 800 //m size bin of the 
PMS particle size distributions in the aircraft measurements (a) on 20 January (lake-effect case) and 
(b) on 22 January (synoptic case). 


3.2.2. Ground-Based Measurement 

[ 42 ] Three ground instruments, laser optical Particle 
Size and Velocity (Parsivel) disdro meter, Geonor weighing 
bucket gauge and 2DVD, were installed at the CARE site. 
The Parsivel disdrometer [Loffler-Mang and Joss, 2000] 
can measure the maximum length of the particle in a hori- 
zontal plane, and its raw output is the number of particles in 
32 size and 32 fall velocity bins. The size range is between 
0.3 to 25 mm, while the falling velocity has a range of 0.05- 


20.8 m s~\ The assumption of measured maximum length 
of the snow aggregate as a diameter results in errors in PSD 
and shape parameters as in the work of Battaglia et al. 
[2010], The Geonor bucket gauge can measure precipita- 
tion rate in the form of liquid accumulation per averaging 
period, typically several minutes. The 2DVD disdrometer 
[Kruger and Krajewski, 2002] measures the maximum width 
and height of falling particle as well as the falling velocity by 
the two orthogonal cameras where the measuring planes 
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have 6 mm apart. This 2DVD system provides the shape 
parameters of snow aggregates to calculate the apparent 
volume and diameter [ Huang et al., 2010], which are used in 
Appendix A. 

[ 43 ] Precipitation PSD is determined to be series of parti- 
cle number concentrations as a function of the maximum 
particle length from the raw output of the Parsivel mea- 
surement. We produce a bulk size-dependent spectrum of 
particle density at the ground level, using the approach that 
is similar to the aircraft measurement analysis. The bulk 
effective particle radius for the Parsivel measurement 
(PCL r e ) is defined following the equation that is similar to 
equation (2): 


32 

X>M 

PCL r e = ^ , (4) 

J= 1 

where r,- is a half of the maximum length in a horizontal 
optical plane of jo, bin and Nj is the particle number con- 
centration of j t ), bin. 

[ 44 ] Solid-phase accumulation volume rate estimated from 
the Parsivel measurement (Ppcl) is combined to the mea- 
surement data of (melted) liquid-phase accumulation mass 
rate (Pgnr) by the Geonor bucket gauge. If it is assumed that 
both Parsivel and Geonor gauge sampled the same accu- 
mulation object, the ratio of the two variables corresponds to 
bulk density (mass per unit volume) of the accumulation. 
The similar approach using 2DVD disdrometer and Geonor 
gauge was employed to determine the relationship between 
density and diameter in the Colorado snowstorm [Brandes 
et al., 2007], The bulk density of surface accumulation, 
GNR-PCL p e , is defined as the ratio of these accumulation 
rates: 


GNR-PCL p e 


Pgnr 

P PCL 


(5) 


These two parameters correspond to representatives of bulk 
particle size and density measured by the ground-based 
instruments. 

[ 45 ] The corresponding PCL r e and GNR-PCL p e are 
calculated from the simulation output using the approaches 
that are similar to those in the aircraft comparison: For all 
surface grids of the output, surface fluxes of 43 binned PSDs 
are calculated using the spatial distributions of PSDs and the 
terminal fall velocities. Then, each 43 binned PSD is dis- 
tributed over the 32 bins identical to the binned size ranges 
of the Parsivel measurement program, and subsequently one 
32 binned PSD is recomposed by totaling. 2DVD correlation 
(Appendix A) is applied to the calculation of the distribu- 
tion. Then, PCL r e is calculated following equation (4). On 
the other hand, GNR-PCL p e is obtained as the ratio of mass 
to volume computed from the 32 binned PSD. The two 
parameters are finally sampled for surface grid points within 
a domain of 4 km square centered at the CARE site for the 
comparison. 

[46] Figures 5b and 5d show scatter diagrams of the 
ground-based observables (PCL r e and GNR-PCL p e ) for 24 h 
(20 and 22 January) from the ground-based measurements 


and the WRF-SBM simulations. Note that the number of 
observed samples (triangle marks) is very small because of 
the limited sampling on the particular spot, in addition to the 
low frequent measurement of the Geonor bucket gauge with 
waiting for melting of solid-phase precipitation; time inte- 
gration of 15 min is employed to ensure a stable measure- 
ment by the bucket gauge in this study. Observed r e -p e 
relationship shows that p e overall decreases with an increase 
in r e , similarly to that in the aircraft measurement. The 
observed plots show smaller effective radii than those of the 
corresponding airborne plots generally, though relatively 
large particles should fall to the surface. The reason for these 
smaller radii is uncertain: probably because of the limited 
sampling, the difference in the definitions of bulk effective 
radii between ground-based and aircraft measurements, time 
integration of sampling, effect of sublimation under cloud 
base, or others. A comparison between the observed plots of 
the two cases suggests that the synoptic case tends to have 
slightly smaller density. This trend is in agreement with that 
of the airborne measurements; that is, the presence or 
absence of riming influences on the microphysical structures 
of surface precipitation also. 

[ 47 ] The corresponding simulation results (Figures 5b 
and 5d) show a part of relatively large density from r e of 
500 pm to 1000 pm in the lake-effect case, whereas there is 
a part of large density around r e of 300 pm in the synoptic 
case. The part of high density in the lake-effect case corre- 
sponds to rimed snow aggregates and the derivative graupel. 
The pattern of the simulated plots is similar to the simulated 
one in the airborne analysis (Figure 5a), except for the large 
size part of r e over 1000 pm. The population of the simu- 
lated plots covers with the observed plots, with the exception 
of the small p e part. A poor assumption of the relationship 
between density and size of nonrimed snow aggregates in 
model microphysics causes the discrepancy at observed 
low density near 0.01 g/cm 3 , as similarly found in the com- 
parison with the aircraft measurement in the synoptic case 
(Figure 5b). 

[48] The ground-based simulated pattern in the synoptic 
case is largely different from those of both the ground-based 
observed one and the airborne simulated one in the same 
case. The part of large density around r e of 300 pm is due to 
a contribution of ice crystals. Although snow aggregates are 
the dominant type of the precipitation, high relative humidity 
for ice over 1 00% and the resultant ice crystals are partially 
simulated near the surface in the vicinity of the CARE site, 
especially after 12Z on the day; the 12Z plot of Figure 2d 
provides a glimpse of the situation. This is considered to 
be mostly due to some forecast errors rather than the prob- 
lem of microphysics. 

3.3. Additional Simulations With Different PBL 
Options or No Riming Process in the Lake-Effect Case 

[ 49 ] The purpose of conducting addition simulations in the 
lake-effect case is to highlight how PBL process affects 
riming process and associated formation of the microphysi- 
cal structure of the snowstorm. In addition, we sought to 
identify a source of the discrepancy between the measure- 
ments and simulation. The specifications of the control and 
addition simulations are summarized in Table 1.(1) Mellar- 
Yamada-Janjic Level 2.5 turbulent closure model (MYJ) is 
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Table 1. List of the Default and Additional Simulations for the 
Lake-Effect Case 


Abbreviation 

Specifications 

MYJ (default) 

Mellar-Yamada-Janjic Level 2.5 


turbulent closure 

MYJ + Nr 

MYJ (No riming 3 ) 

YSU 

Yonsei University PBL b [Hong et al., 2006] 

MRF 

Medium Range Forecast PBL 
model [Hong and Pan, 1996] 


a The riming process is turned off by making riming fraction always zero. 
b PBL, planetary boundary layer. 


used in the run as default. (2) The MYJ PBL model is used 
but the rimed mass of snow aggregates is set always zero, so 
that riming process is turned off and all snow aggregates are 
nonrimed (MYJ + Nr). (3) Yonsei University PBL model 
(YSU) [Hong et ai, 2006]. (4) The Medium Range Forecast 
PBL model (MRF) [Hong and Pan , 1996], 


[so] First, we focus on the difference in the MYJ and 
MYJ + Nr simulations in order to analyze the effect of riming 
process to the characteristics of the lake-effect snowstorm. 
Turn-off of the riming process has largely changed the 
pattern of correlation between particle density and size 
(Figure 8a). The r e -p e relationship in the MYJ + Nr simu- 
lation is roughly logarithmic-linear, and p e uniquely 
decreases with an increase in r e . This pattern is very close to 
that of nonrimed snow aggregates assumed in the SBM 
(Figure la) and simulated in the synoptic case (Figure 5c). 
The part of large density close to p e of 0.4 g/cm 3 has been 
lost because of absence of rimed snow aggregates and 
graupel; this has led to disagreement with the corresponding 
part of the observed scatterplot. However, r e in the 
MYJ + Nr simulation is increased up to excess over the 
maximum value of the observed plots, as compared with 
the result of the MYJ simulation. This increase in particle 
size is considered as being due to replacement of graupels by 
nonrimed snow aggregates with smaller density for the same 
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Figure 8. Same as in the airborne instrument-based diagrams of Figure 5 but for the additional simula- 
tions described in Table 1. Circle marks denote the values derived from the additional simulations. 
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Figure 9. Horizontal distributions of 24 h accumulated surface precipitation from 00:00 UTC 20 January 
to 00:00 UTC 21 January calculated for WRF-SBM (a) Mellar-Yamada-Janjic Level 2.5 turbulent closure 
model (MYJ), (b) MYJ + Nr (no rimming: the riming process is turned off by making riming fraction 
always zero), (c) Yonsei University PBL model (YSU), and (d) Medium Range Forecast PBL model 
(MRF) simulations. 


particle mass. In addition, the minimum value of the simu- 
lated p e approximates to the observed value near 0.01 g/cm 3 . 
This result suggests that the microphysics in the MYJ 
simulation is biased toward overdominance of rimed snow 
aggregates and graupel. Presence of all three types, graupel 
and nonrimed and rimed snow aggregates, is thus necessary to 
represent the widely distributed scatter pattern of r e -p e rela- 
tionship observed by the aircraft measurement (Figure 8a). 
Graupel contributes the part of the observed pattern with 
small r e and large p e \ nonrimed snow aggregates contribute 


the part of the pattern with large r e and small p e , and rimed 
snow aggregates are regarded as the intermediate. 

[ 51 ] The distributions of surface precipitation in the 
MYJ and MYJ + Nr simulations are quite similar (Figures 9a 
and 9b), though the airborne r e -p e relationship has shown a 
remarkable change; the domain-averaged values of the MYJ 
and MY J + Nr simulations are 1 .44 and 1.51 mm in the form 
of the 24 h accumulation, respectively. This result suggests 
that the changes in particle size and density between the 
MYJ and MYJ + Nr simulations has little influence on 
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Lake-effect, WRF-SBM hydrometeor class distribution 


a) YSU 



b) MRF 



Figure 10. Same as Figures 6b and 6d but for the additional simulations described in Table 1. 


determining the transport of hydrometeors to the surface 
through sedimentation in this case. The little change in the 
surface precipitation patterns is due to the small height of the 
simulated snowstorm system and small falling velocities of 
hydrometeor particles. The maximum r ( , of graupel and 
rimed snow in the MYJ run is approximately 700 fim 
(Figure 8a), and the maximum bulk falling velocity of 
graupel is estimated to be approximately 2 m/s (Figure lb). 
In contrast, the corresponding maximum falling velocity of 
nonrimed aggregates in the MYJ + Nr run is estimated to be 
approximately 1 m/s. The maximum difference in horizon- 
tally transported distances is estimated to be approximately 
30 km for descending from 2 km height in the horizontal 
wind with the magnitude of 15 m/s. The actual difference is 
likely smaller than the estimated maximum, because the 
surface precipitation is concentrated at the location of the 
Blue Mountains with the relatively large terrain height. 

[52] Second, we comparatively examined the simulation 
results using different PBL models (MYJ, YSU and MRF). 
The r e of the YSU and MRF simulations are overall increased 
from the MYJ simulation (Figure 8). As a result, the patterns 
of scatterplots of the YSU and MRF simulations correspond 
to be intermediate between those of the MYJ and MYJ + Nr 
simulations. The underestimation of r e in the MYJ run is 
reduced in YSU and especially in MRF run. The pattern of 
r e -p e relationship in the YSU run is more scattered and 
closer to that of the observation, whereas the pattern of the 
MRF run is narrowly distributed and more similar to that of 
the MYJ + Nr run. These r relationships are linked to 
the dominant hydrometer types (Figure 10), as shown in the 
comparative analysis between the lake-effect and synoptic 
cases (Figure 6). The MYJ simulation (Figure 6b) has the 
largest ratios of graupel and liquid droplets among the three 
simulations, whereas the MRF simulation (Figure 10b) has 
the smallest ones. The ratio of graupel is obviously posi- 
tively con-elated with that of liquid droplets, that is, super- 
cooled water. 


[53] The MYJ, YSU and MRF simulations have clearly 
distinct patterns of surface precipitation (Figure 9), as com- 
pared with little difference between the MYJ and MYJ + Nr 
simulations. The domain-averaged values of the YSU and 
MRF simulations are 2.18 and 1.90 mm in the fonn of the 
24 h accumulation, respectively. The YSU and MRF simu- 
lations predict larger amount of surface precipitation, which 
are probably due to the large moisture supply from the lake 
and the difference in the PBL structures (Figure 11). In the 
MRF simulation (Figure 9d), the area of surface precipita- 
tion over 5 mm is extended from the lake Huron to the 
southeast on the downwind side. In addition, the area over 
1 mm hangs over the lake Ontario and extends to outside of 
the plotted domain. On the other hand, in the YSU simula- 
tion (Figure 9c), the area of large precipitation amount over 
12.5 mm is much increased from the MYJ run. The simu- 
lated value at the CARE site is approximated to the observed 
value of 12.3 mm, but the area-wide precipitation seems to 
be much overpredicted as compared with the estimation of 
the surface precipitation from the radar measurement [Shi 
et al., 2010, Figure 3a; Huang et at., 2010, Figure 9]. 

[54] The values of upward moisture flux on the lake Huron 
in the YSU and MRF simulations are increased by an aver- 
age of several tens percentage points from the MYJ run 
(Figure 11). Figure 12 shows comparison of temperature, 
vapor and relative humidity profiles among the simulations 
with the different PBL options at a coastal point of the lake 
Huron. Low temperature and small vapor mixing ratio below 
the altitude of 1.5 km are simulated in the MYJ run, as 
compared with those in the YSU and MRF runs. The PBL 
layer in the MYJ run seems to be thinner, because the 
inversion layer of temperature is simulated at lower level 
than in YSU and MRF runs. These behaviors are considered 
to be due to the difference in upward heat and moisture flux 
on the lake and subsequently related to the difference in 
surface precipitation. However, high relative humidity close 
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Upward moisture flux at the surface (xIO 5 kg/m 2 /s) 
(1200UTC20JAN, Lake-effect case) 





0 0.2 0.4 0.6 0.8 1 1.2 1.4 


Figure 11. Horizontal distributions of upward moisture flux at the surface at 12:00 UTC 20 January 
simulated by the WRF-SBM (a) MYJ, (b) MYJ + Nr, (c) YSU, and (d) MRF runs. 


to 100% is simulated even in the MYJ simulation despite 
smaller vapor mixing ratio, because of lower air temperature. 
Supercooled water is thus generated even in MYJ run in 
spite of the smaller moisture flux from the lake. 

4. Summary and Conclusions 

[ 55 ] In this study, the WRF-ARW 3.1.1 coupled with an 
updated spectral bin microphysics was employed for the 
case study upon cloud microphysical structures of two dis- 
tinct snowfall events (lake-effect and synoptic) over the 
region near the Great Lakes. The updated microphysics 
scheme can predict rimed mass fraction on snow aggregates 
at every grid and time step, so that it can explicitly diagnose 


the bulk density of aggregates and allows a smooth trans- 
formation from snow aggregates to graupel. This imple- 
mentation made it possible to simulate and study detailed ice 
microphysics, considering an influence of riming process. 
The simulation output was compared with the observed data 
by the C-band radar measurement, aircraft probe and 
impactor samplings, and ground-based particle counter 
and bucket gauge samplings. The off-line coupling to the 
G-SDSU allows a direct comparison between the simulation 
and measurement on the primary products base. 

[ 56 ] Our results are summarized as follows: 

[ 57 ] 1. The comparison with the C-band radar measure- 
ment showed that the simulation reproduced the distinct 
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Vertical profiles of T, qv, RH 
(80W 44.5N, 1200UTC20JAN, Lake-effect case) 





Figure 12. Comparisons of vertical profiles of (a) tempera- 
ture (degree C), (b) water vapor mixing ratio (g/kg), and 
(c) relative humidity with respect to water (%) at 44.5°N, 
80°W at 12:00 UTC 20 January, obtained from the WRF- 
SBM simulations with MYJ, YSU, and MRF PBL models. 


basic structures of the two snowfall systems. The lake-effect 
event is characterized by containing localized low and 
intense snowstonn from the lake, whereas the synoptic event 
is characterized by widely distributed high-up and moderate 
snowfall. 

[58] 2. The comparison with the airborne measurements 
showed that the distinct characteristic features of the 


microphysical structures in the two events were successfully 
simulated. The microphysics of the lake-effect snowstorm is 
characterized by rimed snow aggregates and graupel with 
relatively large bulk particle density, while the microphysics 
of the synoptic case snowfall is characterized by nonrimed 
snow aggregates with small particle density. 

[59] 3. Supercooled droplets play a key role in forming 
the characteristic microphysical structure of the lake-effect 
case through riming process. The presence and absence of 
supercooled water in the simulations for both cases were 
justified by the same aircraft measurements. 

[ 60 ] 4. Selection of PBL models has a remarkable impact 
on determining the cloud microphysical structure of the 
lake-effect snowstorm as well as the surface precipitation. In 
contrast, the presence of riming process has a little influence 
on the surface precipitation because of the small height of 
the system. The results of the sensitivity tests also suggest 
that rimed snow aggregates and graupel are overdominant in 
the default run using MYJ PBL model. 

[61] This paper presents a comparative analysis between 
the remotely sensed and in situ measurements and the high- 
resolutional NWP simulation with one of the most sophis- 
ticated microphysical schemes. Detailed discussion on 
effects of riming process on ice cloud microphysics using 
the kind of modeling approach has not been well studied as 
long as we know. This study is an example of great advan- 
tage by introducing the advanced cloud microphysics with a 
gradation of distinct hydrometeor categories to reproduce 
realistic solid-phase microphysical structure in NWP simu- 
lations. The result will be of benefit to the cloud modeling 
community because it may justify the validity of ideal cloud 
resolving simulations using similar advanced microphysics 
previously achieved. Utilizing an up-to-date computational 
server has enabled the high-cost model calculation and the 
off-line handling of the output data including PSDs in great 
volume. 

[ 62 ] In our investigated cases, the two events were clearly 
classified from the view of presence or absence of riming 
process. The condition of the lake-effect case is character- 
ized by the relatively very cold temperature of approxi- 
mately — 20°C at 850 hPa level, which contrasts with the 
corresponding temperature of approximately — 10°C in the 
synoptic case. This cold temperature may be sufficient 
condition for the supersaturation with respect to water 
because of low saturation pressure, but remains to be 
understood how this snowstorm acts in the absence of 
supercooled water at warmer temperature. In addition, it is 
unknown how these distinct snowfall events are influenced 
by difference in the conditions of aerosols acting as con- 
densation nuclei (CN) and/or ice nuclei (IN). The kind of 
clouds grew from a fresh water surface in this lake-effect 
case, whereas similar snowstorms are observed even near 
ocean where sea salt aerosols are enriched in the air mass 
with a strong wind, for example, over the coastal region of 
winter Japan sea [e.g., Tsunogai , 1975]. 

[ 63 ] The change in bulk particle density in the sensitivity 
simulations did not affect the distribution of surface precip- 
itation in the lake-effect case. However, it plays an important 
role to determine the radiative properties that have a large 
influence on the calculation of retrieval algorithms for radar 
and microwave radiometer measurements. For example, 
surface precipitation estimated from a radar measurement 
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may include a large error when the actual particle density is 
different from the value assumed in the algorithm for the 
estimation. Note that this study is involved in the develop- 
ment of the National Aeronautics and Space Administration 
(NASA) synthetic Global Precipitation Measurement (GPM) 
simulator project. A purpose of the project is to offer a set of 
ground validation (GV) constrained 3-D database of WRF- 
SBM output. The database will be used to support testing/ 
developing the prelaunch precipitation retrieval algorithm 
for the measurements by GPM core satellites. The simula- 
tion result of this study has been integrated into the database 
as a preliminary result; therefore, hydrometeor particle size 
and density are important parameters that should be vali- 
dated using observed data. 

[ 64 ] There still remains a large uncertainty in ice micro- 
physics modeling including riming process. Recently, a new 
field campaign was conducted in the same location, so-called 
GPM Cold-Season Precipitation Experiment (GCPEx). It 
deployed more and better in situ and remotely sensed 
instruments to measure characteristics of midlatitude ice 
cloud microphysics. Case studies based on analyses of the 
WRF-SBM simulations and the GCPEx observational data 
sets are subsequently expected. 

Appendix A: The 2DVD Correction Used 
in the Simulator Unit of the Aircraft 
and Ground-Based Measurements Products 

[ 65 ] Snow aggregates are assumed to be spherical shapes 
in the SBM for the reason that it can simplify the modeling 
of microphysical processes such as diffusion growth and 
coagulation. Flowever, maximum dimensions of randomly 
oriented nonspherical aggregates were observed in the in situ 
measurements by aircraft 2-D probes and ground-based 
Parsivel disdrometer. The purpose of the 2DVD-based cor- 
rection is to statistically convert the spherical diameter 
assumed in the SBM into maximum lengths, which corre- 
spond to be the in situ observables, using the 2DVD data. 
Note that this correction is applied to the calculation of the 
aircraft and ground-based measurements in G-SDSU for a 
diagnostic purpose, not within the WRF-SBM weather pre- 
diction simulation. 

[66] The technical description of the 2DVD measurement 
in the C3VP campaign can be found in the work of Huang 
et al. [2010]. The two line scan cameras in the instrument 
provide images of a falling particle from orthogonal hori- 
zontal angles. The shape parameters of particles such as 
shadow areas (Aj and A 2 ), associated vertical maximum 
dimensions (Hi and H 2 ) and horizontal maximum dimen- 
sions (W i and W 2 ), can be obtained from the two snapshots. 
The apparent volume (V app ) and apparent diameter (D app ) 
are defined following equation 1 of Huang et al. [2010] 
using these shape parameters. 

[ 67 ] We sampled the set of 2DVD parameters during 20- 
22 January 2007, and estimated the ratios between maxi- 
mum vertical/horizontal dimensions to the apparent dia- 
meters (R h = H/D app and R w = W/D app ). These ratios are 
normalized for a given D app bins (0-10 mm for each 1 mm 
bin width), as shown in Figure Al. For example, at D app of 
6 mm, horizontal maximum dimension ranges from 4.5 mm 
to 12 mm and the mode value is 7.2 mm. A conditional 
probability density function (PDF) of horizontal maximum 


a) PDF of (vertical diameter)/(apparent diameter) [%] 



b) PDF of (horizontal diameter)/(apparent diameter) [%] 



Dapp [mm] 

Figure Al. Probability density functions of (a) ratio 
between the vertical diameter and apparent diameter and 
(b) ratio between the horizontal diameter and apparent diam- 
eter, from the analysis of the two-dimensional video disd- 
rometer (2DVD) ground measurement. 


dimension for a given apparent diameter can be approxi- 
mated to this PDF in Figure Al. 

[68] If snow aggregates of the SBM are expanded into 
randomly oriented nonspherical particles and the original 
spherical diameter is assumed to be equivalent to D app , the 
PDFs of the maximum vertical/horizontal dimensions are 
approximated to those shown in Figure Al. The maximum 
dimensions are compatible to those actually observed in the 
ground-based Parsivel measurement and approximated in 
the aircraft PMS measurement. The PDF between the max- 
imum (of horizontal and vertical) diameter and the apparent 
diameter is introduced in the calculation of distribution from 
43 binned PSD of the SBM to 33 binned PSD of the simu- 
lated aircraft PMS measurement. The PDF between the 
maximum horizontal diameter and the apparent diameter is 
used in the distribution calculation from 43 binned PSD of 
the SBM to 32 binned precipitation PSD of the simulated 
ground-based Parsivel measurement. This correction is 
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Sphere assumption for snow in place of using 2DVD correction 
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Figure A2. Same as in airborne instrument-based diagrams of Figure 5 but for the simulations without 
the 2DVD correction. 


currently applied only to snow aggregates in the hydrome- 
teor categories of the SBM. 

[69] We briefly discuss the influence of the 2DVD cor- 
rection on the microphysical analysis described in the main 
part of this paper. Figure A2 shows scatter diagrams of the 
airbome-based and ground-based r e -p e relationships without 
the 2DVD correction; the set of the corresponding plots 
using the 2DVD correction is Figure 5. The comparison 
between Figures 5 and A2 shows that the correction has a 
large effect on changing the airbome-based composites, 
whereas little difference is seen in the ground-based com- 
posites. As compared with Figures A2a and A2c, the plots in 
Figures 5a and 5c are biased toward larger r e and remarkably 
smaller p e . This difference is caused by a shift of PSD 


toward the larger size side on the 33 bins by replacing the 
apparent diameter with the maximum diameter in the 2DVD 
correction. The distinctions between the two cases are still 
confirmed, but the agreement with the measurement profiles 
is deteriorated, especially in the airborne composites. This 
result suggests that the assumption of spherical snow body 
in the SBM is not appropriate to be directly compared with 
the airborne or ground-based measurements in this case. The 
difference in degrees of change between the airborne and 
ground-based composites is considered to be due to the 
difference in the PDFs (Figure Al). Note that a large 
uncertainty still remains in this approach, especially in the 
PDF part of the small particle of the diameter less than 
approximately 2 mm because of the instrumental noise. 
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